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Abstract 

The stability of the flow of an incompressible, viscous fluid through a pipe of circu- 
lar cross-section, curved about a central axis is investigated in a weakly nonlinear 
regime. A sinusoidal pressure gradient with zero mean is imposed, acting along the 
pipe. A WKBJ perturbation solution is constructed, taking into account the need 
for an inner solution in the vicinity of the outer bend, which is obtained by identify- 
ing the saddle point of the Taylor number in the complex plane of the cross-sectional 
angle co-ordinate. The equation governing the nonlinear evolution of the leading 
order vortex amplitude is thus determined. The stability analysis of this flow to ax- 
ially periodic disturbances leads to a partial differential system dependent on three 
variables, and since the differential operators in this system are periodic in time, 
Floquet theory may be applied to reduce this system to a coupled infinite system 
of ordinary differential equations, together with homogeneous uncoupled boundary 
conditions. The eigenvalues of this system are calculated numerically to predict 
a critical Taylor number consistent with the analysis of Papageorgiou (1987). A 
discussion of how nonlinear effects alter the linear stability analysis is also given, 
and the nature of the instability determined. 
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1 Introduction 


Our concern is with the stability of the unsteady -viscous flow of an incompressible fluid 
in a pipe of circular cross section, itself curved about a central axis, and subject to a 
sinusoidal pressure gradient of zero-mean. In particular, we are interested in the effects of 
the nonlinear terms of the Navier-Stokes equations on the development of the disturbance 
at the pipe w'alls. The stability of laminar periodic flows, such as in the problem under 
consideration, is often of both mathematical and physical importance. In this case, our 
model is of particular relevance to the fluid mechanics of the cardiovascular system, and 
in particular the possible links to atheroma. Atheroma is increasingly an important 
disease process in middle age, and indeed postmortem results even in young people have 
shown evidence of changes to large arterial vessels. In the initial stages fatty streaks - 
accumulation of lipid in the tunica intima, occur. Subsequently, raised plaques become 
visible on the surface of the arteries. An important feature of atheroma, is that it develops 
at preferred sites in the circulatory system. Early signs may be seen in the major central 
arteries, and not in the peripheral system, and indeed the later onset of atheroma in 
the peripheral circulation might well be due to more advanced lesians in larger vessels 
influencing the blood flow past them and hence the development downstream. The first 
sites, tend to be found near junctions, and curves in the arteries. For example, the 
posterior wall of the descending thoracic aorta almost always has fatty streaks, whilst 
the anterior wall does not. It is thought that one or both of two factors are important in 
the preferential distribution. Firstly the structure of the vessel wall varies so that some 
sites are more prone to development. And secondly that they are subject to different 
influences due to the fluid flow, and this might alter the physiological or biochemical 
processes taking place. Clearly we are interested in the later. Experimental evidence 
indicates that atheroma is linked with low permeability of the endothelium which inhibits 
the efflux of cholesterol. It was initially thought that this low permeability might be due 
to damage to the endothelium due to high wall shear stresses, but in fact it has been 
shown that such damage increase the rate at which molecules may diffuse into and out of 
the plasma. Thus we would expect that less streaks are visible at regions of high shear, 
which agrees with observational evidence. For example, in junctions, streaks are observed 
on the outer walls, but not on the flow divider, and in curved arteries evidence of atheroma 
is more prevalent on the inner bends. For a more detailed account of atherogenesis and 
references to experimental work the reader is referred to Caro, Pedley, Schroter &; Seed . 

Clearly, the model used to study the fluid flow in a curved pipe is far simpler than 
the blood flow within the aortic arch. For example, our model does not take into account 
the distensibility of the arterial wall. However, although this would be of importance in 
determining the local pressure gradient, it does not effect the overall velocity distribution, 
since the wavelength of the pulse of fluid driven by the pressure gradient, is far greater 
than the distance travelled by a typical fluid element in a single cycle. For example, in 
the canine aorta, the pulse wavelength is 3 — 4m, whilst the distance travelled by the 
fluid element is about 10 — 20cm., hence we can assume that during one period, the cross- 
sectional area of the section of the pipe traversed by the fluid is approximately uniform. 
Similarly elastic effects can also be neglected provided that there is no discernible change 
in the pipe cross-section over a length scale of say 10cm. 

Several problems dealing with the stability of oscillatory flows have been investigated, 
both theoretically and experimentally (see Davis (1976)). Periodic laminar flows may be 
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categorised in two ways. Firstly those which have a non-zero mean velocity, in which case 
the disturbance is usually associated with the mean flow and the parameters governing the 
stability problem are dependent on the unperturbed flow field. Secondly, periodic flows 
may be purely oscillatory, that is have zero mean. In general perturbation methods, 
although often applicable to flows of non-zero mean, can not be employed to study 
such disturbances and numerical techniques must be relied upon to analyse the stability 
problem. The flow field investigated in this paper is purely oscillatory in nature and thus 
falls into the latter of these groups. 

Since viscosity is to be included, of importance are the stability mechanisms operating 
within the Stokes layers formed at the pipe walls, and in particular those associated with 
specific types of geometry and local surface behaviour. Since we are investigating the 
flow with a curved pipe, clearly centrifugal effects are of importance, and will effect the 
development of the instability differently depending on whether the surface is locally 
convex or concave. 

The experimental and theoretical investigations of the centrifugal instabilities of a 
Stokes layer in time-periodic flows by Seminara and Hall (1976). They investigated the 
linear stability of the flow induced be a cylinder oscillating harmonically about it’s axis 
in an unbounded fluid. Within the Stokes layer at the surface of the body, the flow is 
shown to be unstable to a Taylor vortex like mode of instability, with the vortices being 
periodic in the azimuthal direction for sufficiently high frequencies of oscillation. Also 
considered was the possible relevance of their study to the stability of the flow within the 
aortic arch. A description of the flow here is complicated, but less important features 
of the problem were neglected. Previous work by Lyne (1971), had shown that, in the 
high frequency limit, viscous effects are confined to a thin layer at the wall, suggesting 
a Stokes layer type flow regime there. Seminara and Hall also investigated the flow 
induced by an oscillating curved pipe, and in particular the stability characteristics of 
the inner (convex) and the outer (concave) bend. They found that at the inner bend the 
disturbance is locally unstable, whilst the outer bend it is locally stable. If, however, the 
flow is driven by an oscillating pressure gradient, as is the case for our problem, the stable 
and unstable regions change positions. A more detailed experimental treatment of this of 
this problem was carried out by Park, Barenghi and Donnelly (1980) and confirmed the 
secondary subharmonic destabilisation of the most dangerous mode found by Seminara 
and Hall , and an approximate description of this breakdown was given by Hall (1981). 

Hall (1984) and Papageorgiou (1987) found that the type of instability mechanism 
investigated be Seminara and Hall can also occur in spatially localised positions in more 
complex boundary layer flows. 

Hall (1984) considered the instability of the two-dimensional flow induced by a trans- 
versely oscillating cylinder, of both elliptical and circular cross-section, in an unbounded 
viscous fluid, in both the linear and weakly nonlinear regimes. For frequencies sufficiently 
high, the cylinder motion drives an unsteady boundary layer which is unstable to Taylor- 
Gortler vortices, localised to regions where the slip velocity of the potential flow, outside 
the boundary layer, is parallel to the motion. For the circular cylinder case in the weakly 
nonlinear regime, it was found that the finite amplitude solutions to the evolution equa- 
tion for the leading order vortex amplitude, bifurcate subcritically from the eigenvalues 
of the linear problem. 

In the problems outlined above, the flow field was essentially two-dimensional. How- 
ever the problem under consideration in this paper comprises of the basic motion in the 
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pipe core, together with a small secondary two-dimensional flow in the cross-section, and 
thus the underlying flow in no longer two-dimensional. The steady problem was first 
studied by Dean (1927,1928), who concluded that the motion depends on a parameter, 
K, defined by, 


where a is the pipe radius, R its radius of curvature about some central axis, and R e the 
Reynolds number. Dean’s analysis was restricted to small values of K, but subsequently 
has been extended numerically to cover moderately large K , and asymptotic boundary 
layer theory has been used to obtain results for very large values of this parameter. 

Until Lyne (1971), the time-dependent problem had not been studied. He gave a 
detailed asymptotic analysis of the fully developed, unperturbed flow in a curved pipe 
under the action of a pressure gradient assumed to be sinusoidal in time, and with zero 
mean. It is also assumed that 6, defined as, 



is small. The flow is found to depend on two parameters, 


W 2 \ 1 _ W 2 a 

auj 2 j H' R$ ~ 


where v is the kinematic viscosity, W is a typical velocity along the pipe, and u> is the 
frequency of oscillation of the basic flow. The parameter t 2 is the ratio of the square of 
the displacement amplitude of the axial motion, (W 2 /au> 2 ), to the radius of curvature, R , 
and is taken to be small. R s may be thought of as the Reynolds number of the secondary 
flow. Also of importance in the analysis is the parameter given by, 


2 _ /2v\ _1_ _ 2<? 
Vo ;) a 2 R* 


(3 is also assumed to be small. Physically, /? represents the ratio of the Stokes layer 
thickness, (2v/u>)?, to the pipe radius, a. Since /3 is small, the viscous effects are confined 
to a thin layer on the pipe wall, whilst the flow in the pipe core is inviscid. The influence 
of the parameter /? on physiological flows was first recognised by Womersley (1955). 

The linear stability problem was investigated theoretically by Papageorgiou (1987) 
and a more detailed discussion of the results obtained is given in §2.3, since this provides 
a basis for the nonlinear problem under consideration here. 

Since the Stokes layer is thin, the streamlines adjacent to the wall can be assumed to 
be of 0(R), and we introduce, T a , the Taylor number for the secondary flow, and defined 
by, 


T = 

1 a — 


w 2 



Rv 2 


2 ~- ( — y 

Rw v \a 2 Lj) 


2 RJ. 


Since we are interested in centrifugal effects, we demand that T a is 0(1), and choose R s 
to be, 

R s = 2 T0~\ 


where T a = 4 T . Note that we may relate the parameters, e, /?, T through the formula, 


e 2 = 0T, 
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and thus the problem may be reduced to one dependent only on the parameters 0 and 
T. Since T a is much smaller than the Taylor number for the axial flow, we may assume 
that the vortices will be aligned with the flow down the pipe and will have characteristic 
wavenumbers based on the Stokes layer thickness, (2v/u)i . In the construction of the 
solution two length scales, 0(1) and 0({3), emerge as being of importance. A possible 
approach, taking into account the different scalings in the pipe cross-section, would be to 
expand the perturbation quantities in powers of /? and apply a WKBJ method. Walton 
(1975) adopted such a method when investigating the narrow spherical annulus problem. 
He anticipated that the critical Taylor number at the equator, above which instabilities 
in the flow may occur, was slightly in excess of that for the corresponding cylinder 
problem. It was found that the solution becomes singular in the vicinity of the of the 
equator, suggesting that an inner expansion is required to smooth out this singularity. 
The solutions to the resulting amplitude equation could be expressed in terms of Airy 
functions, which have the property that, solutions decaying at infinity exhibit oscillatory 
behaviour at minus infinity. Clearly such functions are physically unacceptable, since we 
require that solutions to the velocity distribution must be symmetric about the equator. 
The problem was resolved by Soward and Jones (1983), who identified the Taylor number 
for which inner solutions are well behaved away from the equator. 

For our problem, the linear stability analysis of Papageorgiou showed that, as for the 
narrow spherical annulus problem, the solution to the amplitude equation was singular, in 
this case in the neighbourhood of the outer bend. Again a local inner expansion produced 
physically unacceptable results and the method of Soward and Jones (1983), outlined in 
§2.2, was employed to identify the correct Taylor number for which the inner solution 
is well behaved at ±oo. Papageorgiou ’s investigation also indicated that the instability 
sets in first at the outer bend of the pipe, and it was also shown that there exists no 
critical Taylor number for the corresponding problem at the inner bend, suggesting that 
centrifugal effects are of little importance here. 

The procedure for the remainder of this paper is as follows. In the coming section 
we formulate the problem at hand, obtaining the governing equations for the flow field 
within the Stokes layer , and outline the construction of the inner solution. In §2.3 we 
present a brief summary of the derivation of the linear amplitude equation derived by 
Papageorgiou . This analysis is then extended in §3 to include the nonlinear terms of 
the governing equations and thus the evolution equation for the leading order vortex 
amplitude is obtained in this weakly nonlinear regime. Whilst in §4 numerical results 
for both the linear and nonlinear problems are presented and discussed. Finally we draw 
some conclusions. 


2 Formulation of the Problem 

2.1 Equations of Motion 

Consider the flow of an incompressible viscous fluid in a pipe of circular cross-section 
and radius a, which itself is curved in a circle of radius R about a central axis, as 
illustrated in figure 2.1 The spatial co-ordinates are taken to be, (r, 6,<j>) where r and 0 
are polar co-ordinates within the pipe cross-section, and R<f> measures the distance along 
the pipe. The velocity vector u has components (u, v, w) corresponding to the (r, 0, <j>) 
co-ordinate system and is assumed to be independent of <j), once the basic flow is fully 
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Figure 1: The co-ordinate system 


veloped. 

A sinusoidal pressure gradient of the form, 


d 

<90 



RWu}cos(u>t), 


imposed on the above flow regime. If we consider the Navier-Stokes equations for the 
;tem, and in particular the balance between the viscous work and pressure terms, we 
; that there exists a layer of thickness 0{2vjuj) on the pipe walls inside which the 
cous terms dominate, whilst in the core, we have a potential flow. Similar analysis 
the r- and ^-momentum equations suggests that within the viscous layer, u and v 
: 0(W 2 j3 / Rijj) and 0{W 2 j Ru>) respectively. Thus we introduce the following non- 
nensional notation, 

u ' V „ ,J w 

u ~ w t JrZ'< v ~ wyK,’ w ~ Wi 

p = p(a/R)W* ’ r = o’ r — Loi. 

bstituting (2.1) into the Navier-Stokes equations yields, 

u T + e 2 (uu r + ±vu e - iy 2 ) - w 2 cos 6 = -p T - \P 2 \-§q (u r + yf - , 

v T + e 2 ( uv T + \vv 6 + yuu) + w 2 sin 6 = -\p 9 -(- |/3 2 |: (u r + , 

w T -f e 2 («tu r + hwg) = cos T + \ ft 2 (w rr + pw r + ^w e e) , 

Ur + \u -I- ±Ve = 0. 


( 2 . 1 ) 

(2.2a) 

(2.2b) 

(2.2c) 

(2. 2d) 


: ease in notation, we have ignored the superscripts on the new, non-dimensional 
iables. The system of equations above, (2.2a-d), describes the viscous flow field both 
ir the wall of the pipe and in its core. The continuity equation of (2. 2d) may be 
isfied by introducing a non-dimensional stream function 0, with u and v given by, 

u = -ipe, v = -0 r . 

r 

; are primarily interested in the stability of the flow field adjacent to the walls, i.e. 
,hin the Stokes layer, which is of thickness, (2v/u))5 /3. Thus we introduce the following 
lings for r and 0, 

t, 1 = - r), = /? -1 0, 
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where r, and T are the new radial co-ordinate and stream function respectively, inside 
the Stokes layer. The solution for the basic flow within this layer, is given by Lyne, as a 
series expansion in /?, 

T = $0 + /?Ti + /3 2 ^2 + • • • , (2.3a) 

w = wbo + 0 w Bi + I^ 2 wb2 + • • • • (2.3b) 

In (2.3a,b), T, and WBi for i = 0, 1,2. . ., are functions of t,t, and 0 and are found for 

fixed values of R s , the secondary Reynolds number. Expressions for To, Ti may be found 

in Lyne (1971). As was stated in the introduction, we require the asymptotic solution, 
as ft — » 0 and R s — ► oo. 

The first two terms in the expansion for wb are given by, 


wbo = sin r-e ,, sin(r-7;), (2.4a) 

wbi = ~\ r ! e ~ 1] sin(r — r ,) . (2.4b) 

Denoting the basic flow by, Ug = (ug,ug,u>g,pg), wo have the following expressions 
for ub and vg, 


wg — d u Bo + /3 2 ubi + - • • , (2.5a) 

vb = vbo + Pvbi + - . • , (2.5b) 


where, 



VBi = 


dV t 

dr, ’ 


( 2 . 6 ) 


for i = 0, 1 

The solutions, (2.5a,b), are strictly valid only in the limit as R s — > oo with (3 held 
fixed. However, since R s has no explicit effect on the equations for T 0 , wbo etc within the 
Stokes layer, but only effects the core flow, where it acts like the conventional Reynolds 
number, the solutions above, (2.4a, b), (2.5a, b) are indeed valid representations of the 
flow field as (3 — *■ 0. 

Before proceeding with the analysis, we re-write the Navier- Stokes, equations (2.2a-d), 
in terms of the Stokes variable, r, and thus obtain the system of equations governing the 
flow field within the Stokes layer. 


2.2 Construction of the Inner Solution 

In the introduction we mentioned that, for both the narrow spherical annulus near the 
equator and the curved pipe problem in the vicinity of the outer bend, inner expansions 
lead to physically unrealistic results. Such problems can be resolved by the analytic 
continuation of the solution into the complex 0-plane. Soward and Jones first used such 
a method to resolve the narrow spherical annulus problem, and Papageorgiou modified 
their method when consideri ng the linear stability of the flow in a curved pipe. We now 
describe this approach when applied to the study undertaken here. 

We wish to obtain a dispersion relation for the Taylor number, T of the form, 

T = T(k,0,cr), (2.7) 
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where a is the growth rate, and k, 9 are as previously defined. T is an analytic function 
of the complex variables k , 0, <x, but we require that, for physical flows, it remain real and 
constant for all 0. At a minimum of T, for physically acceptable solutions, the following 
conditions must hold, 

T k = 0, T e = 0. (2.8) 

In Hall (1984), these conditions are satisfied on the real 0-axis, but clearly complex 
values of 9 are also permissible. For the narrow spherical annulus problem, studied by 
Walton, Soward and Jones found that, at the minimum value of T on 9 = 0, it was found 
that although Tk = 0, T\ did not vanish, and the resulting amplitude equation is of the 
Airy-type, which does not lead to physically realistic solutions. 

If we are to obtain valid solutions that describe the vortex amplitudes at the critical 
Taylor number, we must locate the saddle points of T at which the conditions of (2.8) 
are met. Since, for our problem, this does not occur for real 9 we must analytically 
continue the solution into the complex 0-plane. Suppose that such a point exists, say 
at 0 = 0o, then the solution found in the neighbourhood of this point provides both the 
inner solution, that removes the singularity at the outer bend, 0 = 0, and also, a set of 
asymptotic solutions that match with those valid away from the saddle point, and which 
describe the flow field as /? — * 0, for all real 0. 

Due to the symmetry of the problem, we shall be considering only modes that are 
symmetric with respect to 0. Thus, 

Re (0 O ) = 0, 7m (fc 0 ) = 0. (2.9) 

We can obtain two linearly independent asymptotic solutions to the problem, of a 
WKBJ type which take the form, 

q± = exp {±ip~ l j fc(0')d0'l qi ± (0,?7,r) + 0{/3). (2.10) 


o(tf ) 



p. s (9=0) p + 


Figure 2: A diagram illustrating the domains into which the complex 0-plane is divided. 
The inner solution is valid within the shaded area, a circle of radius 0(/3 2). 

The stability problem, under consideration, is an eigenvalue problem and has a solu- 
tion only if a condition of the following form, 

F{k,9,T) = 0, (2.11) 
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is satisfied. The dependence of F on the variables 0, T follows from the search for 
neutrally stable solutions, that is solutions with zero growth rate, and by a Fourier 
analysis of the solutions in time. The condition that T* = 0 of (2.8), implies that there 
exists a repeated root of (2.11) at 9 = 0 O . Let this root occur when, 

k { 1 ) ( 6 0 ) = kW( 6 0 ) = k 0 . (2.12) 


Thus, the point P, at which 0 = 0 O in the complex 0-plane defines a transition point, in 
the neighbourhood of which the WKBJ solutions, (2.10), become identical. The second 
condition Te = 0 implies that, in the limit as 0 — *■ 0 the transition points coincide, at 
0 = 0o- Hence there exists a pair of anti-Stokes lines dividing the complex plane into four 
regions. These are defined by, 


Im {k w (9') - k { 1 ) ( 0 ')) d9 ' | = 0. 

Note that we may also define the corresponding Stokes lines by, 

Re jjf* (f- l2 W') - M ' | = 0. 


(2.13) 


(2.14) 


Figure 2-2, is a schematic representation of how the complex 0-plane is divided by the 
Stokes and anti-Stokes lines. The exact contours of the principle curves have not been 
calculated from (2.13). A more detailed description is given by Soward and Jones (1983). 
The reader is also referred to Drazin and Reid (1981). 

Since we are concerned with only the symmetric modes, the Stokes line, PS is the 
imaginary 0-axis, whilst the anti-Stokes lines, PP- and PP+ join 0 = 0 O to 0 = —9\ and 
0 = +0i respectively. Hence we can divide the complex 0-plane into the regions Ri — ► R 4 . 

We define qi and q 2 to be the asymptotic solutions corresponding to the roots j feW 
and kW respectively. Let us assume, that q x be the dominant solution in region R\. That 
is the error associated with this solution, is 0(/?qi), which is small compared with bfq x 
with in this region. On crossing the anti-Stokes line PP L into R 2 , the associated error 
changes to 0 ( 0 q 2 ), which is nolonger small compared with the solution q 1? and qi is 
said to recessive in i? 2 . In i?3, q! remains recessive, and on crossing the anti-Stokes line 
PP+ into i?4, qi again provides a good approximation to the exact solution. Similarly, 
q 2 is dominant in regions R 2 + R 3 and recessive in R\ + R 4 . For a full description of 
phase-integral methods the reader is referred to Heading (1962). 

It is, in general, possible to approximate the true solution uniformly along the real 
0-axis by, 

qi 0 > 0 

qi — ( q 2 0 < 0 ' ' 

where ( is a Stokes constant, found by obtaining an inner solution valid for | 0 — 0 O 
O (/?^) and matching this with the outer solution. 

WKBJ solutions corresponding to roots of (2.11) other than and k ^ are not 
required, since they are not linked to the transition point at 0 = 0 O . We require that the 
Stokes constant, ( vanishes, in order that qt is the approximate solution along the whole 
of the real 0-axis. At S, 6 — 0, the wavenumber and its derivative, with respect to 0, 
have the following properties, 




Im(k s ) = 0, Im(k S 9 ) > 0. 


(2.16) 
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Thus the disturbance takes its maximum value at S. The asymptotic solutions, (2.10) 
are valid provided that | 9 - 0 o |» (3 * , and control the matching with the inner solutions. 


2,3 Review of the Linear Problem 

We now present a brief summary of the inner analysis and the results obtained by Papa- 
georgiou (1987). 

An inner variable 0 is defined by, 

6-6 0 = fiiQ. (2.17) 


and note that, 0 is of 0(1). 

We also introduce the following notation for the perturbation vector, q 


„ dv dw 
drj ’ dr] 


91= ,—,u,v,w 


The perturbation vector, q is then expanded in powers of /?2 as, 

q = d!(0)qi£ + d 2 (Q)q 2 E + d 3 (0)q 3 £/? + . . . + c.c., 


(2.18) 


(2.19) 


where, E = j ik 0 Q//3^ and c.c. denotes the complex conjugate. The q,-, * = 1,2,... are 
defined in a similar manner to q. Having introduced a disturbance, q, into the basic flow 
within the Stokes layer, we describe the total flow as, 


{u,v,w,p) = ( u B ,v B ,w B ,p B ) + ei(u,i),u;,p), 

where is a small vanishing parameter, and the subscript, B is, once again, used to 
denote quantities of the unperturbed, basic flow field. 

It is now a straight forward procedure to substitute the total flow field into the 
governing equations and linearise with respect to ei. The leading order problem for qi 


reduces to, 

Co(qi) — o, 

(2.20) 

where Co is 

the operator defined by, 



< 

1 

<! 



C° — I— Ao + B— . 

or] dr 

(2.21) 


In the above, Ao and B are 6x6 matrices not given here. 
At 0(0^), q 2 satisfies, 


= cb©Miqi + OdjlVbqi + diM 3 qi. (2.22) 

Here M t , i = 1,2, 3 are 6x6 matrices not given here. The solution to (2.22) is of the form, 

d 2 q[2 = -idieq^ + iQd iq^ 2) + (2.23) 


T <9q 2 i . , Tjdcu 

V ' Aoq2 + b D7 
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where a^i = 1,2, 3 are evaluated using the following system of equations, 

Co(q2^) = — C'A; 0 o(qi), 

Co(q2 2) ) = iC© o 0 (qi), 

C 0 (q 2 3) ) = *'Cfcoo(qie). (2.24) 

In (2.24), the subscripts, k 0 , 9 0 denote partial differentiation with respect to k 0 , 0 o respec- 
tively, and zero represents evaluation at 6 = 0 O . 

The evolution equation for the leading order amplitude function d\ is determined by 
imposing a solvability condition on the differential system obtained at 0(0). The required 
condition is, 


+ I/v0di© -(- lpd\ + IgdiQ + l/?0di + l50 2 di = 0, (2.25) 

where the coefficients, Tv — > I5 are double integrals over p- and r- space, for example, 

Im = jH £* o V T (M (0) qx + drdr,. ( 2 . 26 ) 

3 Formulation of the Nonlinear Stability Problem 

We now go on to describe how nonlinear effects alter the linear stability analysis outlined 
in §2.3. The terms fundamental, mean and first-harmonic refer to the dependence of 
the disturbance on the ^-co-ordinate. We recall that the 9- dependent amplitude of 
the disturbance satisfies (2.25) and that exponentially decaying solutions exist for only 
certain values of the Taylor number, T, that is those satisfying the conditions (2.8). We 
scale the amplitude of the disturbance in such a way as to introduce non-linear effects 
which modify the linear amplitude equation. In order to retain the linear structure of 
(2.25), we expand the Taylor number as, 

T = T oc + 0T u (3.1) 

where T 0 c , is the critical value of T 0 . We first assume that the velocity field perturbation 
is O(0 S ). At O(0 2S ), the fundamental interacts with itself, through the nonlinear terms of 
the Navier-Stokes equations, to produce first-harmonic and mean flow correction terms. 
These then interact with the fundamental, reinforcing it at O(0 3S ). Thus for (2.25) to be 
modified by nonlinear effects we must choose 6 = 

We now outline the construction of the inner solution following the method of Papa- 
georgiou (1987), described briefly in §2.3 . We define an 0(1), inner variable 0 by, 

6 - 0 O = 0/3* , (3.2) 

where 0 = 0 O is the value of 0 , for which the conditions, (2.8) are satisfied. k 0 and T 0 are 
similarly defined. The disturbance quantities, u,v,w and p are then expanded as, 

u = Ednitn cos <f>0* + [£Jd 2 i«2i cos <t> + £ 2 d 22 U22 cos 2 <f> + d 20 w 2 0] 0 
+ [^d 3 iu 3 i cos 4> + E 2 d 32 u 32 cos 2 <f> + E 3 d 33 u 33 cos 3 <j> + d 30 u 30 ] /?* + ... + c.c., (3.3a) 
v = Ed n v u cos 4>0* + \Ed 2 \v 2 \ cos <j> + E 2 d 22 v 22 cos 2<p + ^20^20] 0 
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+ \Ed 3 iVz\ cos 0 + E 2 d 32 v 32 cos 20 + E 3 d 33 v 33 cos 30 + d 30 U3o] /?? + ... + c.c., (3.3b) 
w = Ed n w u cos 00% + \Ed 2 \w 2 i cos <j> + E 2 d 22 w 22 cos 20 + ^20^20] 0 
+ ]Ed 3 iw 3 1 cos 0 + E 2 d 32 w 32 cos 20 + E 3 d 33 iv 33 cos 3 0 + ^30^30] 0* + ... + c.c.(3.3c) 
p = Ednpu cos 00% + [^c?2iP2i cos <f> + E 2 d 22 p 22 cos 2 4> + c? 2 oP2o] 0 
+ [Ed 3 X p 3l cos 0 + E 2 d 32 p 32 cos 20 + E 3 d 33 p 33 cos 30 + c? 3 oP3o] /?* + ••• + c.c., (3.3d) 

The coefficients in the above expansions depends only on the variables, r, 0 and 77 , and 
c.c. denotes the complex conjugate . E is given by, 


E = exp 



Hence, the total flow may be expressed as, 


(3.4) 


(u,V,W,p) = ( U B ,VBiW B ,PB ) + («,t7,t&,p), 


(3.5) 


where the subscript B denotes the basic flow. At this point we introduce the perturbation 
vector q tJ , defined as, 

r dv ‘> 1 (3.6) 


q ij = 


Pij ? 


, Uij , Vij 7 


’ < 9 t / ’ dp 

We proceed by substituting, (3.5) into the governing equations and successively equate 
like powers of 02 , so that at 0 ( 0 2) we find that the fundamental, qn, satisfies the linear 
stability problem, with 0, k and T evaluated at the saddle point of T, 


Co(qn) — 0, 


where, C 0 is the operator defined by, 


r\ r\ 

C 0 = I— - Ao + B— . 
Op Or 


(3.7) 


(3.8) 


Ao is the matrix A evaluated at the saddle point of T and B is the 6x6 matrix given 
below, 


A = 


0 —\ik 0 \k 2 + ikT 0 VBo 


2 ik 

0 

0 

0 

0 


2 

0 0 
0 0 
0 0 


0 


—2wbq cos 9 


-2Tqv B ot, k 2 + 2ikT 0 v B o 4te B0 sin 0 
-2T 0 w B o n 0 k 2 + 2ikT 0 v B Q 


1 

0 


0 

1 


0 

0 

0 


ik 

0 

0 


0 

0 

0 


(3.9) 


B 


0 0 0 -1 

0 0 0 0 

0 0 0 

0 0 0 

0 0 0 

0 0 0 


0 

0 

0 

0 


0 

-2 

0 

0 

0 

0 


0 

0 

-2 

0 

0 

0 


(3.10) 
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At 0(0), the fundermental, q 2 i, satisfies, 

d 2 iCo (q 2 i) = duMiqn© -f dii@M 2 qn + dn©Miqn (3-11) 

where the matrices M„ i = 1,2 are related to .the operator C 0 through the following 
equations, 


A* -dCo 

Ml = 'dk 0 

(.3.12a) 

s 

to 

II 

1 

Oil 2? 

£>19' 

(3.12b) 

And we may obtain a solution to (3.11), by writing, 


<^21^21 = -«d 11 ©qj> 1 i ) + iQd n c$ + rf n q^, 

(3.13) 

where q 2 i , i = 1, 2, 3 are to be determined. Substitution of (3.13) into (3.11) and equating 

coefficients of d n ©, 0dn and dn yields the following equations for 

q^? and q 2 3) , 

Co(q 2 ^) = — C* 0 o(qn), 

(3.14a) 

Co(q2i^) = *C© 0 o(qn), 

(3.14b) 

Co(q 2 i ) ) = *C*; o o(qxi0). 

(3.14c) 


In (3.14a, b,c), the subscripts k o ,0 0 denote partial differentiation with respect to k 0 ,9 0 
respectively, and the zero subscript denotes evaluation at the saddle point of T, that is 
where, 9 = 9$. In addition to the fundamental mode, first harmonic, q 22 , and mean flow 
correction, q 20 terms are also generated. After some manipulation, we find that these 
maybe expressed as, 


d 22 q22 — ^n422) (3.15a) 

^20*120 — diidnq 2 o. (3.15b) 

In addition, q 2 2, q 2 o satisfy, 

C 0 (q 2 2) = D, (3.16) 

Co (q 22 ) = D. (3-17) 

In the above equations for q 22 and q 2 o, D and D represent the terms arising form the 
interaction of the leading order fundamental with itself. 

The amplitude function, d n is obtained be imposing a solvability condition on the 
differential system for the fundamental, q 31 at 0 ( 0 2), 

dz\ C 0 (q 3 i) = due© [M<°)q„ + M^V] 

+dn©0 [N (0) q n + N (1) q5V + N (2) q< 2) ] 

+d n [P (0) qn + P (2) q^ + P (3) q^] 

+dn@ [Q (0) qn + Q (1) q 2 V + Q (3) q^] 

+e<*„ [R (0) qn + R< 2 >q< 2) + R< 3 )qg>] 

+0 2 d n [S<°>q„ + S( 2 )q( 2) ] 

+d n | d 2 x | T^ 0) qn, (3.18) 
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where, M (,) , N^, P^, Q^, R^, S^, i — 0,1,2, 3 are 6x6 matrices not given here, and 
T (0) , contains the terms arising from the interaction of the leading order fundamental 
with the first harmonic and mean flow correction terms. By considering the partial 
differential system adjoint to (3.7), we find that. -a solution to the above equation exists 

if, 

dn©elM + Odn©I/v + diilp + dnelg + 0dnlp + 0 2 dnls + d\\ | d\ x | Ip = 0, (3.19) 


where I M etc are double integrals over r/- and r-space. For example, 

Im = jH /^V T [M (0) q n + M (1) 4V] dr dr/. 

And V satisfies the adjoint differential system, 

Cj(V) = 0 

Vi = v 2 = V 3 = o at T) = 0 


where Cj is the operator, 


**-£ +a * t+bt !- 


(3.20) 

(3.21) 

(3.22) 

(3.23) 


The amplitude equation (3.19) must be solved subject to the boundary conditions, 

dn —> 0 as | 0 |— *• oo. (3.24) 


4 Numerical Solution and Results 


In this section, we describe the numerical procedure used to solve the leading order 
eigenvalue problem subject to boundary conditions of no-slip at the pipe wall and decay 
of the disturbance at infinity. 

The leading order problem must be solved subject to conditions on k 0 and 6 0 . The 
eigenrelation is solved for fixed values of k and 0 until the value of T 0 corresponding to 
physically acceptable solutions is located. 

On the basis of Floquet theory, since the basic flow is periodic in time, we assume 
that the disturbances are also periodic and carry out a Fourier expansion in time for the 
perturbation quantities. For neutrally stable solutions this takes the form, 

OO 

q=GEqV" T + c.c, (4.1) 

— OO 


where G is a constant and q n are functions of rj alone. After some manipulation, the 
equations for the leading order problem may be reduced to a pair of coupled partial 
differential equations for u and to, 


d>_ 

dr] 2 


-k 2 




kj Un + 2ikTo (k 2 v B o + itBOriTj) 


«n 2tkT\)V sou lirfri 


-Ak 2 WBoW U cos 0 o - 4 ik sin 0 o (wB 0 r,Wn + ^Bo^iir,) = 0, 


d 2 

dr ] 2 


- k 2 



+ 2ToWboti u ii ~ 2ikT 0 VBoWn — 0 - 


(4.2a) 

(4.2b) 
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Substitution of (4.1) into the above system would lead to an infinite set of coupled 
ordinary differential equations at successive powers of e' nT . In order to reduce this to a 
more tractable system of equations, we set u n = w n = 0 for | n |> M. In addition, we 
replace oo by where 77 ^ and M are chosen to give the degree of accuracy required. 
For large 77 , (4.2a-b) reduces to, 

(jhf ~ ^ 2 ) Un ~ \*kTok 2 sm0 o u u + ^iTTosinfloUiir,*, 

—4 k 2 sin Tien cos 0 O — 4 ik sin 6 0 sin tw u v = 0, 
w u + +^ikT 0 sin = 0. 

The above system, has three independent solutions with the correct behaviour as 77 — > 00 . 
These may be used to integrate (4.2a-b) from 77 = 7700 to 77 = 0, in this instance using a 
fourth order Runge-Kutta scheme with h steps in the interval of integration, to provide 
3(2M + 1) solutions. Combining the independent solutions of the reduced system at 77 = 0 
we may satisfy all but one of the boundary conditions, and the remaining condition is 
met if T(k,0 ) is an eigenvalue of the system. 




(4.3a) 

(4.3b) 



Figure 3: The eigenfunctions tqj 1 , u n 3 and plotted as a function of 77 for the critical 
value of 7 q 


Our calculations, for M = 6 , 77 ^ = 12.5 and h = 251 gave results of T 0 = 10.7301, 

9 0 — 0.3525 i and k 0 = 0.5313, comparing favourably with those given by Papageorgiou (1987). 
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The eigenfunctions corresponding to (4.2a-b) were found to have the property that 

u"j = 0 n even, (4.4a) 

tt>jj =0 n odd. (4.4b) 

The functions, ujj 1 and evaluated at the critical values of k 0 ,6 0 and To are illustrated 
in figure 3, together with the eigenfunctions for wj" 3 - The solutions were normalised such 
that the first derivative of was set to 1.0. 

The solution to the adjoint solution was calculated in the same manner and used as a 
check on the values found above, since the eigenvalues of the adjoint system are identical 
with those of (4.2a-b). The functions, Vj" 1 and V® of the adjoint system are illustrated 
in figure 4. The adjoint eigenfunctions are such that, 



Figure 4: The eigenfunctions V x 1 and V° for the adjoint problem corresponding to the 
critical Taylor number, To 


V" = 0 n even, (4.5a) 

V" = 0 n odd. (4.5b) 

In addition solutions for the fundamental, mean flow correction and first harmonic 
at 0((3) were calculated in a similar fashion and hence the integral coefficients of the 
amplitude equations (2.25,3.19) could be calculated, in this case using Simpson’s rule. 

Before proceeding we first simplify the linear amplitude equation, so that the solution 
may be expressed in terms of known functions. Using the transformations, 

d x = e( _V(C “ 7)2 )e(- A7) Z/, (4.6a) 
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C = 0 + 7, 


(4.6b) 


we may re-write (2.25) as 


where 


Z" + p 


H + Cp T Ti - 




-C)z, = o, 


n = 

H = 
Cpt = 


A 2 -7V + ^ 

lM 
IPT 
Im ’ 


a 



/m ) 



(4.7) 


(4.8a) 

(4.8b) 

(4.8c) 


In the above transformations, a, A and 7 are suitably chosen constants, depending on the 
integral coefficients of (3.19). Since we have assumed in (4.1) that the coefficients of the 
Fourier expansion depend only on r), we find that A and 7 and both zero. The amplitude 
function d iy and hence, Zi, must decay as, | 0 |— » 00 , and hence the solutions to ( 4 . 7 ) 
are, 

(0 = Z* (C) = U n (- (n + 1) ,2^ic) , (4.9) 

where U n , is the n th parabolic cylinder function, and the corresponding value of Ti is, 


T\ — Tin — 


2jj} (n + \)-H 


Cpr 


(4.10) 


The functions Zi n (() have n — 1 zeros for ( € (— 00 , + 00 ) and depending on the odd/even 
nature of n, are odd or even in (. All the functions tend to zero like exp(-fi *C 2 /2) and 
the least stable mode corresponds to the n = 0 case, when 


Z i0 = exp 



5 


and 


Ti = Tic = 


fA -_H 

Cpr 


(4.11) 


Our calculations showed that Tic = 3.005. 

Returning to the nonlinear amplitude equation (3.19), using the transformations, 
(4.6a-b), this maybe re-written in the following form, 


where 


z" + f, _ f 2 j 2 = _ CtZ \z\*F((), (4.12) 

F (C) = | | a , (4.13a) 

cv = (4.13b) 

16 




Figure 5: The numerically calculated bifurcating solution of (4.12) 


and i/, //, Cpt are given by (4.8a-c). 

If we expand Z and Tj in terms of some small parameter e, 

Z = £2 Z\ + e* Z 3 + ... , 

T x = Tc + tTp + . . . , 

and substitute these into (4.12) then we find at leading order, Z\ satisfies the linear 
amplitude equation (4.7) and hence, 


Z\ — <f>Zi, 

where 4 1 is an arbitrary constant determined at higher order. In addition, Tc = Tic , 
where Tie is given be (4.11). At 0(e*), we find that Z 3 satisfies, 

Z” + H ( — + - c 2 ) = -C^Zi I <f>Zi I 2 F (0 - CptTp^Zx (4.14) 

In order for solutions to the above equation to exist, a solvability condition must be 
satisfied. By this we mean that the right hand side of the differential equation must be 
orthogonal to the solution of the adjoint problem to (4.7). It should be noted that (4.7) 
is self-adjoint. The required condition is, 

(4 - 15) 

C-T J 2 
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where J,, i = 1, 2 are given below, 


/ OO 

I Zi | 2 ac, 

-oo 

/ OO 

l zi i 4 ^ «)<*<■ 

-oo 

If we assume that | T\ — 7c |<C 1 then (4.15) maybe be written in the form, 

1 z |2 ~ ,Tl _ Tc) 1 z ' |2 (4,16) 



Figure 6: Comparison of solutions for T x = 2.9. The solid line represents the solution the 
full problem, (4.12) and the o symbol the approximate solution, (4.14) 


Since Cpt/Ct was found to be positive, finite amplitude solutions of (4.12) only 
exist locally near T\ = Tic for T\ < T\q. Thus the solutions to (4.12) will bifurcate 
subcritically from the eigenvalues of the linear problem. The subcritical nature of the 
bifurcation was confirmed by numerically integrating (4.12) using a shooting procedure. 
We note that, since the differential operator in (4.12) is even in and F (() is also 
an even function for A = 0 thus we can expect that the solution Z(() is either even 
or odd in £, depending on the conditions applied, with even solutions corresponding to 
Z'( 0) = 0 and odd solutions to, Z'( 0 + ) = Z'( 0~), and in this case the former condition 
was applied. The results are shown in figure 5, where we have plotted the amplitude of 
the first mode, evaluated at ( = 0, as a function of 7\. It is possible that higher-order 
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Figure 7: Comparison of | Z( 0) | 2 from (4.12) - solid line, (4.14) - small dashed line and 
(4.16) - dashed line 


nonlinear effects may reverse this result, producing supercritical equilibrium solutions, 
however higher-order calculations would be required in order to substantiate/disprove 
this conjecture. 

Subcritical solutions are known to be unstable, and by allowing the amplitude function 
Z to be dependent on some slow time variable this may be shown. Thus the fully nonlinear 
problem must be solved numerically in order to find out how the flow develops. 

In figure 6 we compare the results obtained from numerically integrating (4.12) - solid 
line, with those from the expansion procedure carried out in the neighbourhood of Tic ■ 
We took, T\ = 2.9, and as illustrated the results were found to be almost identical. 

The results obtained by numerically integrating the nonlinear amplitude equation, 
(4.12) - solid line, are compared with those obtained from the expansion about Tic, - 
small dashed line, and the approximate value for | Z | 2 from (4.16) - dashed line, in figure 
7. As expected the results from the two approximate methods are in good agreement 
with those from the numerical integration in the neighbourhood of Tic- 


5 Conclusions 

We have obtained the equation governing the nonlinear evolution of the leading-order 
vortex amplitude function dn(0). As expected, the linear terms of this equation are of 
the same form as those of the equivalent equation found from the linear stability analysis. 
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Indeed, we would also expect that the coefficients of the linear terms, I I N , Ip, Ig, Ip, Is, 
will be identical to those for the linear evolution equation, found by Papageorgiou (1987), 
taking account of the differences in notation. 

The results of our numerical calculations predict that the instability is subcritical 
in nature, and thus we presume that close to the critical Taylor number, sufficiently 
large perturbations to the basic state will grow. These disturbances might tend to some 
equilibrium, due to stabilisation by higher-order nonlinear effects, and thus some form 
of steady state reached. Alternatively, even if no stabilisation takes place, due to the 
localised nature of the instability, we might expect that some periodicity along the pipe 
exists. For example, Tollmien-Schlichting waves are subcritical in nature, but can be 
observed in parallel or nearly parallel flows. The fully nonlinear problem would need to 
be solved numerically in order to investigate the subcritical nature of the bifurcation, 
and find the flow into which the disturbance evolves. 

In the introduction, we mentioned the relevance of the model to the study of the fluid 
mechanics of blood flow in large arteries and in particular the aortic arch. Typically for 
the canine cardiovascular system, the ascending aorta is 1.5cm in diameter, the mean 
velocity is approximately, 20 cms _1 ,/? s ~ 4000, ~ 0.1, and 8 has a value of about 0.2. 

Our analysis is not inconsistent with these values, assuming that the mechanics of the 
blood flow are not significantly altered in the limit 8 — > 0. 

However, one important feature of physiological flows is that in general the pressure 
gradient has a non-zero mean component, and thus the Dean number, D must be included 
as a parameter in the problem. For the canine aorta this has a value of « 2000. The 
problem under consideration here, becomes physiologically viable if we consider the case, 
(3-+0 ,D < R s , (see Papageorgiou ). In this case, the flow-field within the Stokes layer, 
to leading order, is described by Lyne’s analysis and the stability of the solution is as 
described here. 
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